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SUMMARY: Finite element (FE) analysis was performed on 3 -point and 4-point 
bending test configurations of ninety degree oriented glass-epoxy and graphite-epoxy 
composite beams to identify deviations from beam theory predictions. Both linear and 
geometric non-linear analyses were performed using the ABAQUS ? finite element code. 
The 3 -point and 4-point bending specimens were first modeled with two-dimensional 
elements. Three-dimensional finite element models were then performed for selected 4- 
point bending configurations to study the stress distribution across the width of the 
specimens and compare the results to the stresses computed from two-dimensional plane- 
strain and plane-stress analyses and the stresses from beam theory. Stresses for all 
configurations were analyzed at load levels corresponding to the measured transverse 
tensile strength of the materials. 

For 3-point bend test configurations, both the linear and geometric non-linear 2D plane- 
strain and plane-stress analyses yielded similar results. The maximum tensile stresses 
under the center load nose calculated from the FE analysis were slightly lower than 
stresses predicted by beam theory. The difference (maximum of 4%) was greatest for the 
shortest span analyzed. 

For 4-point bend test configurations, both the plane-stress and plane-strain 2D linear 
analysis results agreed closely with beam theory except right below the load points. 
However, 2D geometric non-linear analyses deviated slightly from beam theory 
throughout the inner span as well as below the load points. Plane-stress results deviated 
from beam theory more than plane-strain results. The maximum tensile stresses between 
the inner span load points were slightly greater than the beam theory result. This 
difference was greatest (maximum of 4%) for configurations with the shortest spans 
between inner and outer load points. A contact analysis was also performed in order to 
investigate the influence of modeling the roller versus modeling the support as a simple 
boundary condition at one nodal point. A configuration with the shortest span between 
inner and outer load points was modeled for the 24-ply and 36-ply IM7/8552 layups and 
the 24-ply S2/8552 layup. Generally, for all configurations investigated, the discrepancy 
between the FE and beam theory results became smaller (max. 2%) when the rollers were 
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modeled in conjunction with contact analysis. Hence, the beam theory yields a reasonably 
accurate value for the maximum tensile stress in bending compared to 2D FE analysis. 

The 3D linear FE analysis of the 4-point configurations agreed closely with beam theory, 
except right below the load points. The 3D linear FE results at the specimen edge agreed 
with 2D plane-stress results, and the 3D linear FE results in the center of the specimen 
agreed with 2D plane-strain results. The 3D geometric non-linear analyses deviated 
slightly from beam theory throughout the inner span as well as under load points. The 3D 
geometric non-linear FE results at the specimen edge agreed with the 2D plane-stress 
results. For the 12.7 mm (0.50 in.) wide IM7/8552 specimens, the 3D geometric non- 
linear FE results in the center of the specimen agreed with 2D plane-strain results. 
However, for the 6.35 mm (0.25 in.) wide S2/8552 specimens, the 3D geometric non- 
linear FE results in the center were less than 2D plane-strain results, indicating that these 
specimens were not wide enough to achieve full constraint. 

The utility of the FE results is primarily for guidance in the choice of beam thickness, 
width, and configuration. For the 3-point bend configuration, longer spans are preferred 
to minimize the error in beam theory data reduction. Similarly, for the 4-point bend 
configurations, a longer span between the inner and outer load noses, at least equal to the 
span between the inner load noses, results in less error compared to beam theory. In 
addition, these FE results indicate that the span between the inner load noses should not 
be too long to avoid obtaining a non-uniform maximum stress between the inner load 
noses. Finally, the 3D analysis indicates that specimens should be sufficiently wide to 
achieve a fully constrained state of plane-strain at the center of the specimen width. 


KEYWORDS: transverse tensile strength, matrix cracking, finite element analysis, glass 
epoxy, graphite epoxy, bending tests 

INTRODUCTION 

Matrix ply cracking is a common initial damage mechanism in fiber reinforced 
composites. Matrix ply cracking alone is seldom catastrophic for laminates subjected to 
membrane loading. However, for composite structures that undergo bending, or other 
out-of-plane loading, the formation of matrix cracks may lead to immediate catastrophic 
delamination formation and growth [1,2], Hence, accurate characterization of the 
transverse tensile strength of composite materials is needed to accurately predict matrix 
ply cracking in these structures. 

Previously, 90 degree unidirectional glass-epoxy laminates [3] and carbon-epoxy 
laminates [4] were tested in three and four point bending to characterize composite 
transverse tensile strength. The influence of edge flaws due to machining, and surface 
flaws due to manufacture and handling, were assessed by testing specimens in their as- 
manufactured and machined condition. In addition, specimens with polished edges, 
and/or bottom failure surfaces, were tested. Subsequent 3-point and 4-point bend tests of 
glass-epoxy lamina and carbon-epoxy lamina [5] were also performed under cyclic 
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loading. In each of these studies, the maximum tensile stress at failure under monotonic 
and cyclic loading was determined using beam theory expressions. 

The goal of this investigation was to study under which conditions simple beam theory 
yielded sufficiently accurate results when used to calculate the stresses in the 3-point and 
4-point bending specimens. For the entire investigation, the ABAQUS 8 finite element 
software was used. For comparison, both linear and geometric nonlinear analysis 
procedures were used. The three-point and four-point bending specimens were first 
modeled entirely with two-dimensional elements. Then, three-dimensional finite element 
models were used to study the stress distribution across the width of the specimens and 
compare the results to the stresses computed from two-dimensional plane-strain and 
plane-stress analyses and the stresses from simple beam theory. 

Materials 

FE analyses were performed on two materials, S2/8552 glass-epoxy [3,5] and IM7/8552 
carbon-epoxy [4,5], Material properties used in the analysis are given in Table 1. The 
S2/8552 specimens were cut from a 24-ply panel. The average panel thickness was 5.56 
mm, corresponding to an average ply thickness of 0.231 mm. The IM7/8552 specimens 
were cut from both 24-ply and 36-ply panels. The average panel thickness was 3.29 mm 
for the 24-ply panel and 4.94 mm for the 36-ply panel, corresponding to an average ply 
thickness of 0.137 mm. 


Specimen Configurations 

Three-point bending tests were evaluated in three configurations (A, B, and C) 
corresponding to span lengths of 25.4, 50.8, and 76.2 mm (figure 1(a)). Four point 
bending tests were evaluated using the four configurations (A 1 , A2, B2, B3) shown in 
figure 1(b). 


Beam Theory 


For the three point bending configurations shown in figure 1(a), with span length, s, and 
width, b, the tensile bending stress distribution at the bottom of the beam is shown in 
figure2a. Specimen strengths were determined in references 3&4 from the beam theory 
expression for the maximum tension stress, G max , under the center load nose using the 
maximum load at failure, P c , and the average specimen thickness, ta Vg , 


a 


max 


3P c s 

2bt - g 


(i) 


For the four point bending tests shown in figure 1(b), with outer span, s, inner span, s-f, 

and width, b, the tensile bending stress distribution at the bottom of the beam is shown in 
figure 2b. Specimen strengths were determined in references 3&4 from the beam theory 
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expression for the maximum uniform tension stress, o max , between the inner load points 
using the maximum load at failure, P c , and the average specimen thickness, t avg , 


a 


max 


3P c l 

2bt avg 


( 2 ) 


Finite Element (FE) Analysis Formulation 

Two typical two-dimensional finite element models of the specimens, boundary 
conditions and loads applied in the simulations are shown in Figure 3. A typical three- 
dimensional finite element model of the specimens used in the simulations is shown in 
Figure 4. For both the 2D and 3D analyses, the ABAQUS* input file was generated 
automatically using a FORTRAN routine (see Appendix 1). This user- written routine 
required only the input of the geometrical data, the material property data given in Table 
I, and the mesh parameters shown in Figure 3. The use of this routine allowed the quick 
generation of a large number of models for the different specimen configurations and 
simplified the otherwise tedious mesh generation and refinement. 

The two-dimensional specimen cross sections were modeled using eight-noded 
quadrilateral plane-strain (CPE8R) or plane-stress (CRS8R) elements using quadratic 
shape functions and a reduced (2x2) integration scheme. These elements were selected to 
avoid shear locking in bending, which is common in first-order, fully integrated elements, 
such as CPE4 and CPS4 [6], The numerical formulation of first order elements gives rise 
to shear strains that do not really exist. Therefore, these elements are too stiff in bending, 
and many elements over the thickness are required to obtain acceptable results. Elements 
where a lower-order, reduced integration are used to form the element stiffness, such as 
the CPE8R and CPS8R elements used in this study, usually require fewer elements 
through the thickness, provide more accurate results in bending, and yield reduced 
running times. Loads were prescribed at single nodal points in the models. However, for 
configuration B2 that exhibited the greatest deviation from beam theory, rollers were also 
included in the model at the outside load points (see Appendix 2). 

For the 3D analyses, specimens were modeled with ABAQUS solid, twenty-noded, 
C3D20R hexahedral elements using quadratic shape functions and a reduced (2x2x3) 
integration scheme. The C3D20R element utilizes reduced integration to form the 
element stiffness, and hence, provides accurate results in bending which yield reduced 
running times compared to other 3D elements [6,7], 

Computation of Stresses 

Averaged stresses at nodes were computed along the bottom (tension) surface of the 
beam specimens where the highest tensile stresses occur. The results were extracted for 
additional post-processing directly from the ABAQUS s binary result file using a user 
written FORTRAN routine (See Appendix 1 ) [6,7], 
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Two Dimensional FE Analysis Results 
Three Point Bending 

Figure 5 shows the finite element mesh for the S2/8552 glass epoxy material, 3-point 
bend configuration A (figure la). The mesh is shown in the deformed position at the 
failure load with the stress contours superimposed illustrating the compression stress near 
the point of load application at the top of the beam and the high tensile stresses at the 
bottom of the beam under the center load nose. Tension stress results at the bottom of the 
beam will be quantified in later figures. Similar meshes were generated for the other 3- 
point bend configurations in figure la using the technique described earlier. 

The open symbols in figure 6 show the stresses calculated from the 2D, plane-strain, 
linear finite element analysis, Ofe, normalized by the maximum tension stress calculated 
from the beam theory, Obt ( O max in equation 1). For comparison, the normalized stress 
distribution from beam theory is superimposed on the finite element results. This 
comparison indicates that the peak tensile stress calculated from the finite element 
analysis is slightly (4%) less than the beam theory prediction. Similar results were 
obtained for the 2D plane-stress linear finite element analysis, as well as the 2D 
geometric non-linear plane-strain and plane-stress FE analyses. In order to establish the 
accuracy of this comparison, the 2D geometric non-linear analysis was also performed 
using a refined mesh as shown in figure 7. This analysis yielded more stress results in the 
critical area under the load nose. However, this refined mesh yielded essentially no 
change in the maximum stress calculated directly under the load nose, which was also 
less than the beam theory (figure 8). 

Figure 9 summarizes the peak stress comparisons for the 24-ply S2/8552 glass-epoxy 3- 
point bend configurations in terms of percent A a, where 

Ao = aBT ~ aFE xlOO (3) 

Obt 

As shown in figure 9, the greatest difference between all the 2D FE analyses and the 
beam theory was about 4% for the shortest span, configuration A. The longer span 
configurations had difference between 1 % and 2%. 

Figure 10 summarizes the peak stress comparisons for the 24-ply IM7/8552 
carbon-epoxy 3-point bend configurations. The greatest difference between all the 2D FE 
analyses and the beam theory was about 2% for the shortest span, configuration A. The 
longer span configurations had differences of 1 % or less. Figure 1 1 summarizes the peak 
stress comparisons for the 36-ply IM7/8552 carbon-epoxy 3-point bend configurations. 
The greatest difference between all the 2D FE analyses and the beam theory was about 
3% for the shortest span, configuration A. The longer span configurations had difference 
between 1 % and 2%. 
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The results from figures 10-11 indicate that the beam theory results become less accurate 
for thicker beams of a given span, and for shorter spans of the same beam thickness. This 
is the same trend expected for the significance of transverse shearing deformation in 3- 
point bending of beams. This shearing deformation is not reflected in the beam theory 
equations 1 and 2. 

Four Point Bending 

Figure 12 shows the finite element mesh for the S2/8552 glass-epoxy material, 4-point 
bend configuration A1 (figure lb). The mesh is shown in the deformed position at the 
failure load with the stress contours superimposed illustrating the compression stress 
between the two inner load noses at the top of the beam and the high tensile stresses at 
the bottom of the beam under the inner load noses. Tension stress results at the bottom of 
the beam will be quantified in later figures. Similar meshes were generated for the other 
4-point bend configurations in figure lb using the technique described earlier. 

The open symbols in figure 13 show the stresses calculated from the 2D plane-strain, 
linear finite element analysis, Gfe, normalized by the maximum tension stress calculated 
from the beam theory, Gbt ( O max in equation 2 ). For comparison, the normalized stress 
distribution from beam theory is superimposed on the finite element results. This 
comparison indicates that the peak tensile stress calculated from the finite element 
analysis is slightly less than the beam theory prediction just below the load nose as 
observed in the 3-point bend case. Flowever, the 4-point bend FE results become slightly 
greater than the beam theory stress just inside the load nose before assuming the beam 
theory result in the majority of the span between the inner load noses. Similar results 
were obtained for the 2D plane-stress, linear finite element analysis. The 2D plane-strain 
geometric non-linear FE analysis yielded a similar FE stress distribution, except the 
plateau in the majority of the span between the inner load noses was larger (1.5%) than 
the beam theory result ( figure 14). The 2D plane-stress, geometric non-linear FE analysis, 
yielded a similar FE stress distribution, but with an even greater discrepancy (2.5%) from 
the beam theory between the two center load noses (figure 15). In order to establish the 
accuracy of this comparison, the 2D geometric non-linear analysis was also performed 
using a refined mesh as shown in figure 16. This yielded more stress results between the 
inner load noses, but there was essential no change in the plateau of the distribution 
between the inner load noses (figure 17). Figures 18 a-d compare the 2D geometric non- 
linear plane-strain and plane-stress results in terms of percent Ag (equation 3) between 
the inner load noses for the four S2/8552 glass-epoxy 4-point bend configurations 
(figurelb). For all four configurations, the plane-stress results differ from the beam 
theory predictions by 0.5 to 1.5% more than the plane-strain results. 

The 2D plane-stress geometric non-linear FE results in terms of percent Ag (equation 3) 
between the inner load noses are summarized for all four S2/8552 glass epoxy 4-point 
bend configurations in figure 19. For the plateau in the majority of the span between the 
inner load noses, the greatest discrepancy between the FE and beam theory results occurs 
for configurations B2 (4%) and A1 (2.5%). These two configurations have the shortest 
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span between the inner and outer load noses (12.7 mm). Figure 19 also illustrates the 
non-uniform stress distribution between the inner load noses for the two configurations, 
B2 and B3, with the greatest inner span, s-l = 50.8 mm. These same trends were noted 

for the 24-ply (figure 20) and 36-ply (figure 21) IM7/8552 carbon-epoxy, 4-point bend 
configurations. The thinner 24-ply laminate had slightly lower discrepancies (only 3% for 
configuration B2, for example) than the thicker 36-ply laminate of the same material. 
However, the non-uniform stress distribution between the inner load noses was more 
pronounced in the thinner laminates. 

Contact Analysis 

A contact analysis was performed in order to investigate the influence of modeling the 
roller versus modeling the support as a simple boundary condition at one nodal point (see 
appendix 2). It was expected that accounting for the roller would change the contact 
location and result in a shorter span, s, of the specimen and thus result in a different stress 
distribution along the bottom surface of the specimen. The roller support was modeled 
and corresponding contact analysis was performed for specimen configuration B2, only, 
because discrepancy in stresses between the finite element and the beam theory results 
were most pronounced for this configuration as shown in Figures 19-21. The B2 
configuration was modeled for the 24-ply and 36-ply 1M7/8552 layups and the 24-ply 
S2/8552 layup. As shown in the appendix, for all configurations investigated the 
discrepancy between the FE and beam theory results became smaller (max. 3%) when the 
rollers were modeled in conjunction with contact analysis instead of simple supported 
boundary conditions. Hence, the beam theory yields a reasonably accurate value for the 
maximum tensile stress in bending compared to 2D FE analysis. 

Three Dimensional FE Analysis results 

The S2/8552 glass-epoxy material, 4-point bend configuration B2 (figure lb) was chosen 
for the 3D analysis since it exhibited the greatest deviation from beam theory in the 2D 
non-linear analysis. Figure 22 shows the finite element mesh for a 6.35mm wide S2/8552 
glass-epoxy, 4-point bend configuration B2. The mesh is shown in the deformed position 
at the failure load. The 3D model had the same in-plane mesh refinement as the 2D 
analyses, but with ten elements across the specimen width. 

Figure 23 shows the 3D geometric non-linear FE stress results in terms of percent Ag 
( equation3) between the inner load noses and across the specimen width. Near the load 
noses, the stresses peak near the center of the specimen width. However, between the 
load noses the stress peaks near the edges of the specimen. Ideally, the stresses 
calculated near the specimen edges should be in good agreement with the 2D plane-stress 
results and the stresses calculated near the center of the specimen should be in good 
agreement with the 2D plane-strain results. Figure 24a compares the 3D stresses at the 
specimen edge to the 2D plane-stress results for both the linear and geometric non-linear 
cases. Good agreement was found in both cases, with only the non-linear case showing 
the significant deviation from the beam theory result between the inner load noses. Figure 
24b compares the 3D stresses at the center of the specimen width to the 2D plane-strain 
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results for both the linear and geometric non-linear cases. Good agreement was found for 
the linear case, with no significant deviation from the beam theory result between the 
inner load noses. However, the 3D non-linear case indicated a greater deviation from the 
beam theory result between the inner load noses (4%) than did the 2D plane-strain result 
(3%). This indicated that these beams were not quite wide enough to achieve a fully 
constrained (plane-strain) condition at the specimen center. 

To investigate the influence of specimen width further, a 3D analysis was also performed 
for a 12.7mm wide, 24-ply, IM7/8552 carbon-epoxy, 4-point bend configuration A1 
(figure lb). This configuration was used in reference 5 to characterize transverse tension 
fatigue life. Figure 25 shows the finite element mesh for the 24-ply, IM7/8552 graphite- 
epoxy material, 4-point bend configuration Al. The mesh is shown in the deformed 
position at the failure load. The 3D model had the same in-plane mesh refinement as the 
2D analyses, but with ten elements across the specimen width. 

Figure 26 shows the 3D geometric non-linear FE stress results in terms of percent Ac 
(equation 3 ) between the inner load noses and across the specimen width. Figure 27a 
compares the 3D stresses at the specimen edge to the 2D plane-stress results for both the 
linear and geometric non-linear cases. Good agreement was found in both cases, with 
only the non-linear case showing a slight deviation (1.5 %) from the beam theory result 
between the inner load noses. Figure 27b compares the 3D stresses at the center of the 
specimen width to the 2D plane-strain results for both the linear and geometric non-linear 
cases. Good agreement was found for the linear case, with no significant deviation from 
the beam theory result between the inner load noses. Furthermore, for this 12.7 mm wide 
specimen, there was also good agreement between the 2D non-linear plane-strain solution 
and the 3D non-linear solution at the center of the specimen width, indicating that a fully 
constrained (plane-strain) condition had been achieved at the specimen center. Both the 
2D and 3D non-linear analyses results indicated a slight (1.5%) deviation from the beam 
theory at the specimen center. 

Analysis of Results 

Figures 28-30 compare the 2D stress analysis results for 3-point bend configurations A, 
B, and C for the S2/8552, 24-ply IM7/8552, and 36-ply IM7/8552 materials, respectively, 
to the histogram of failure locations measured for these materials [4], The histograms 
show the location and frequency of failures for the number of replicates tested of a given 
configuration. These figures indicate that most failures occur near the center of the span, 
but not always directly below the central load nose where the tensile stress is a maximum 
(figure 2a). As noted in references 3-5, all failures do not occur directly below the load 
nose because the failure is sensitive to the presence of flaws in the specimen 
microstructure. Both the beam theory and the finite element analysis assume that the 
material is uniform and homogeneous throughout. Hence, it would be meaningless to 
reduce the data based on the FE stress at a given span location relative to the center load 
nose because the failure is due to the local stresses associated with micro-structural flaws. 
Rather, in order to characterize the strength and the associated volume dependence, data 
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are fit to either a normal or Weibull distribution ranked according to the calculated 
maximum stress at failure from the beam theory [4,5]. 

Figures 31 a-d compare the 2D stress analysis results for 4-point bend configurations Al, 
A2, B2, and B3, respectively, for the S2/8552 glass epoxy material to the histogram of 
failure locations measured for these materials [3,4], These figures indicate that most 
failures occur near the center of the span, between the inner two load noses where the 
tensile stress is a maximum (figure 2b). The remaining failures occur below the inner 
load noses as was observed in the 3-point bend specimens. 

Figures 32 a-c compare the 2D stress analysis results for 4-point bend configurations Al, 
A2, and B2, respectively, for the 24-ply IM7/8552 carbon-epoxy material to the 
histogram of failure locations measured for these materials [4], Figures 33 a-d compare 
the 2D stress analysis results for 4-point bend configurations Al, A2, B2, and B3, 
respectively, for the 36-ply IM7/8552 carbon-epoxy material to the histogram of failure 
locations measured for these materials [4], Unlike the glass-epoxy material, a greater 
percentage of failures occur below the inner load noses for all the configurations. In 
reference 4, this was attributed to the large variability in specimen thickness observed in 
the graphite material, which is not accounted for in either the beam theory or FE analysis, 
but leads to a biased loading towards one, or the other, inner load noses in the 
experiments. 

The utility of the FE results is primarily for guidance in the choice of beam thickness, 
width, and configuration. For the 3-point bend configuration, longer spans are preferred 
to minimize the error in beam theory data reduction. Similarly, for the 4-point bend 
configurations, a longer span between the inner and outer load noses, at least equal to the 
span between the inner load noses, results in less error compared to beam theory. In 
addition, these FE results indicate that the span between the inner load noses should not 
be too long to avoid obtaining a non-uniform maximum stress between the inner load 
noses. Finally, the 3D analysis indicates that specimens should be sufficiently wide to 
achieve a fully constrained state of plane-strain at the center of the specimen width. 

Conclusions 

For 3-point bend test configurations, both the linear and geometric non-linear 2D plane- 
strain and plane-stress analyses yielded similar results. The maximum tensile stresses 
under the center load nose calculated from the finite element (FE) analysis were slightly 
lower than stresses predicted by beam theory. The difference (maximum of 4%) was 
greatest for the shortest span analyzed. 

For 4-point bend test configurations, both the plane-stress and plane-strain 2D linear 
analysis results agreed closely with beam theory except right below the load points. 
However, 2D geometric non-linear analyses deviated slightly from beam theory 
throughout the inner span as well as below the load points. Plane-stress results deviated 
from beam theory more than plane-strain results. The maximum tensile stresses between 
the inner span load points were slightly greater than the beam theory result. This 
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difference was greatest (maximum of 4%) for configurations with the shortest spans 
between inner and outer load points. A contact analysis was also performed in order to 
investigate the influence of modeling the roller versus modeling the support as a simple 
boundary condition at one nodal point. The B2 configuration was modeled for the 24-ply 
and 36-ply IM7/8552 layups and the 24-ply S2/8552 layup. Generally, for all 
configurations investigated, the discrepancy between the FE and beam theory results 
became smaller (max. 2%) when the rollers were modeled in conjunction with contact 
analysis. Hence, the beam theory yields a reasonably accurate value for the maximum 
tensile stress in bending compared to 2D FE analysis. 

The 3D linear FE analysis of the 4-point configurations agreed closely with beam theory, 
except right below the load points. The 3D linear FE results at the specimen edge agreed 
with 2D plane-stress results, and the 3D linear FE results in the center of the specimen 
agreed with 2D plane-strain results. The 3D geometric non-linear analyses deviated 
slightly from beam theory throughout the inner span as well as under load points. The 3D 
geometric non-linear FE results at the specimen edge agreed with the 2D plane-stress 
results. For the 12.7 mm (0.50 in.) wide IM7/8552 specimens, the 3D geometric non- 
linear FE results in the center of the specimen agreed with 2D plane-strain results. 
However, for the 6.35 mm (0.25 in.) wide S2/8552 specimens, the 3D geometric non- 
linear FE results in the center were less than 2D plane-strain results, indicating that these 
specimens were not wide enough to achieve full constraint. 

The utility of the FE results is primarily for guidance in the choice of beam thickness, 
width, and configuration. For the 3-point bend configuration, longer spans are preferred 
to minimize the error in beam theory data reduction. Similarly, for the 4-point bend 
configurations, a longer span between the inner and outer load noses, at least equal to the 
span between the inner load noses, results in less error compared to beam theory. In 
addition, these FE results indicate that the span between the inner load noses should not 
be too long to avoid obtaining a non-uniform maximum stress between the inner load 
noses. Finally, the 3D analysis indicates that specimens should be sufficiently wide to 
achieve a fully constrained state of plane-strain at the center of the specimen width. 
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Appendix 1 

ft® 

A general overview of the analysis procedure is given in Figure Al. The ABAQUS 
input file was generated automatically using a user written FORTRAN routine 
(genbeam.f in figure Al), which, in addition to the material data given in Table I, only 
required the input of the loading, the geometrical data and mesh parameters as shown in 
Figures 3 and 4. The use of this routine (flow chart shown in Figure A2) allowed the 
quick generation of a large number of models for the different specimens and simplified 
the otherwise tedious mesh refinement and convergence studies. 

The stresses along the bottom side/surface of the beam specimens were of prime interest 
as this is the location where the highest tensile stresses are to be expected. Therefore, the 
averaged stresses at nodes were extracted for additional post-processing directly from the 
ABAQUS ' binary result file using another user written FORTRAN routine (stress. f in 
figure Al) the flow chart of which is depicted in Figure A3. 
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Appendix 2 

Finite Element (FE) Analysis Formulation 

A contact analysis was performed in order to investigate the influence of modeling the 
roller versus modeling the support as a simple boundary condition at one nodal point. It 
was expected that accounting for the roller would change the contact location and result 
in a shorter span, s, of the specimen and thus result in a different stress distribution along 
the bottom surface of the specimen. The contact analysis was performed using only two- 
dimensional finite element models of the specimens as shown in Figure A4a. Only 20% 
of the outer part of the roller was modeled with finite elements as shown in detail in 
Figure A4b. A roller with a diameter of 4.73 mm was modeled which corresponds to the 
roller diameter of the test fixture used in references 3&4. Material properties for the steel 
rollers used in the simulation are given in Table 1. Displacements in x and y direction 
were suppressed for all nodes along the inner radius of the model and along the line of 
symmetry of the roller to simulate the rigid behavior of the roller. Contact was modeled 
between two deformable bodies with small relative motion without friction using the 
concept of master and slave surfaces available in ABAQUS ® [6]. The top surface of the 
elements in the refined section of the roller modeled made up the master, or target , 
surface as shown in Figure A4b, while the surface of the elements in the refined section 
at the bottom of the specimen defined the slave or contact surface used in the analysis. 

Contact Analysis results 

The roller support was modeled and corresponding contact analysis was performed for 
specimen configuration B2, only, because discrepancy in stresses between the finite 
element and the beam theory results were most pronounced for this configuration as 
shown in Figures 19-21. Additionally, the thinner 24-ply IM7/8552 B2 configuration 
exhibited a more pronounced non-uniform stress distribution between the inner load 
noses than the thick laminates. Therefore, the B2 configuration was modeled for the 24- 
ply and 36-ply IM7/8552 layups and the 24-ply S2/8552 layup. 

Figure A5 shows the finite element mesh for the 24-ply IM7/8552 carbon/ epoxy material 
B2 type specimen. The full specimen and the details are shown in the deformed position 
at the failure load. As shown in the detail of Figure A5, contact occurs only locally. Due 
to the local rotation of the specimen at the support, the contact location is changed and 
moves inward which reduces the effective span s by about 0.6 mm at each support. 

The tension stresses computed at the bottom of the specimen, using 2D contact analysis 
where the roller support were modeled, are shown in Figures A6 a-c. The difference 
between the beam theory and FE stresses for all configurations are lower compared to 
results obtained from analysis where simple supported boundary conditions had been 
used. This may be explained by considering equation (2). Because the distance between 
the load points (s-f) remains unchanged in the simulation, the reduction of the effective 

span, s, causes the same reduction in the length, l , which is used in equation (2) to 
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calculate the stresses. For all configurations investigated, the discrepancy between the FE 
and beam theory results became smaller (max. 3%) when the rollers were modeled in 
conjunction with contact analysis instead of simple supported boundary conditions. 


TABLE I. MATERIAL PROPERTIES. 


S2/8552 Unidirectional Glass-Epoxy Prepreg 

An = 47.71 GPa 

A 22 = 12.27 GPa 

A 33 = 12.27 GPa 

V 12 = 0.278 

V 13 = 0.278 

V 23 = 0.403 

Gn - 4.83 GPa 

Go = 4.83 GPa 

G 23 = 4.48 GPa 

IM7/8552 Unidirectional Carbon-Epoxy Prepreg 

An = 161.0 GPa 

A 22 = 11.38 GPa 

A 33 = 11.38 GPa 

v 12 = 0.32 

vn = 0.32 

V 23 = 0.436 

Gi 2 = 5.17 GPa 

Gi 3 = 5.17 GPa 

G 23 = 3.98 GPa 

Steel 

A = 210.0 GPa 


v = 0.3 
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(a) 3 -point bending (b) 4-point bending 

Figure 1. Bending test configurations. 


► i/2 \* -s-i ►) i/2 \* 

P P/2 P/2 

■ ■ 



(a) 3 -point bending (b) 4-point bending 

Figure 2. Beam theory tension stress distribution on the bottom of beam. 
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t: specimen thickness 

nr: number of elements in refined section 

nsp number of elements in overhang 

ns2: number of elements in section between support and load point 
nt: number of elements over the thickness 

(a) Two dimensional FE model of three point bending specimen 




s-L: span between inner load points 

ns3: number of elements in the center section 

(b) Two dimensional FE model of four point bending specimen 


Figure 3. Two dimensional FE models and corresponding load and boundary conditions 
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L: total specimen length 
s: span 

r: length of refined section around 
supports and load points 
t: specimen thickness 
b: width 



nr: number of elements in refined section 
nsi: number of elements in overhang 
ns2: number of elements in section 
between support and load point 
nt: number of elements over the thickness 
nb: number of elements across the width 
s-L: span between inner load points 
ns3: number of elements in the center section 


Figure 4. Three dimensional FE model of four point bending specimen. 
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17 


Fig. 5 Finite Element Mesh showing Deformed Geometry and Bending Stresses 
for S2/8552 Glass Epoxy 3-Point Bend Configuration A 
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Fig. 6 Normalized Tension Stress Distribution from 2D Plane-Strain Linear FE Analysis 

S2/8552 Glass-epoxy 3-point bending Configuration A 
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Fig. 7 Finite Element Fine Mesh showing Deformed Geometry and Bending Stresses 
for S2/8552 Glass Epoxy 3-Point Bend Configuration A 


Beam Theory 
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Fig. 8 Normalized Tension Stress Distribution from 2D Plane-Stress non-linear FE Analysis 
Fine Mesh for S2/8552 Glass-epoxy 3-Point Bending Configuration A 
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Fig. 9 Comparison of Discrepancy between Beam Theory and FE analysis 
S2/8552 Glass-epoxy, 3-point bending configurations 
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Fig. 10 Comparison of Discrepancy between Beam Theory and FE analysis 
24-ply IM7/8552 Carbon-epoxy, 3-point bending configurations 
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Fig.l 1 Comparison of Discrepancy between Beam Theory and FE analysis 
36-ply IM7/8552 Carbon-epoxy, 3-point bending configurations 
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Fig. 12 Finite Element Mesh showing Defonned Geometry and Bending Stresses 
for S2/8552 Glass Epoxy 4-Point Bend configuration A1 



FE Analysis 
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Fig. 13 Normalized Tension Stress Distribution from 2D Plane-Strain Linear FE Analysis 

S2/8552 Glass-epoxy 4-point bending Configuration A1 



FE Analysis 
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Fig. 14 Normalized Tension Stress Distribution from 2D Plane-Strain non-linear FE Analysis 

S2/8552 Glass-epoxy 4-point bending Configuration A1 



FE Analysis 
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Fig. 15 Normalized Tension Stress Distribution from 2D Plane-Stress non-linear FE Analysis 

S2/8552 Glass-epoxy 4-point bending Configuration A1 
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Fig. 16 Finite Element Fine Mesh showing Deformed Geometry and Bending Stresses 
for S2/8552 Glass Epoxy, 4-Point Bend configuration A1 


FE Analysis 
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Fig. 17 Normalized Tension Stress Distribution from 2D Plane-Stress non-linear FE Analysis 
Fine Mesh for S2/8552 Glass-epoxy 4-point bending Configuration A1 
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(b) Configuration A2 
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(c) Configuration B2 (d) Configuration B3 

Fig. 18 Comparison of Discrepancy between Beam Theory and Plane-Stress and 
Plane-Strain FE Analysis for S2/8552 Glass-epoxy 4-point bending configurations 
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Fig. 19 Comparison of Discrepancy between Beam Theory and 2D Plane-Stress 
non-linear FE analysis for S2/8552 Glass-epoxy 4-point bending Configurations 
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Fig. 20 Comparison of Discrepancy between Beam Theory and 2D Plane-Stress 
non-linear FE analysis for 24-ply IM7/8552 Carbon-epoxy 4-point bending Configurations 
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Fig. 21 Comparison of Discrepancy between Beam Theory and 2D Plane-Stress 
non-linear FE analysis for 36-ply IM7/8552 Carbon-epoxy 4-point bending Configurations 
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W) 

E 



(a) 2D Plane Stress vs. 3D Edge results 



(b) 2D Plane Strain vs. 3D Center results 


Fig. 24 Comparison of Discrepancy between Beam Theory and 2D and 3D 
FE Analysis for S2/8552 Glass-epoxy 4-point bending Configuration B2 
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Fig. 26 Beam Theory vs. 3D Finite Element - IM7/8552 Configuration A1 




Length, mm 


(a) 2D Plane Stress vs. 3D Edge results 



Length, mm 


(b) 2D Plane Strain vs. 3D Center results 

Fig.27 Comparison of Discrepancy between Beam Theory and 2D and 3D 
FE Analysis, 24-ply IM7/8552 Carbon-epoxy 4-point bending Configuration A1 
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Fig .28 Comparison of Normalized Bending Stress and Failure Location 
S 2/855 2 3-point bending Configurations 
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Fig .29 Comparison of Normalized Bending Stress and Failure Location 
24-ply IM7/8552 3-point bending Configurations 


41 



Conf. A 


Conf. B 



-10 -5 0 5 10 


Distance from center load point (mm) 



Distance from center load point (mm) 



Fig .30 Comparison of Normalized Bending Stress and Failure Location 
36-ply IM7/8552 3-point bending Configurations 
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Fig .31 Comparison of Bending Stress Distributions and Failure Locations, S2/8552 4-point bending 
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Fig. 32 Comparison of Bending Stress Distributions and Failure Locations, 24-ply IM7/8552, 4-point bending 
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Fig. 33 Comparison of Bending Stress Distributions and Failure Locations, 36-ply IM7/8552, 4-point bending 



/ user specified \ 

/ input file \ 

/ genbeam.data \ 

\ contains loads, / 

\ dimensions and / 
\ material data / 


genbeam.f 

read input data 
calculate nodal point 
coordinates, 
position of loads and 
boundaries 

write ABAQUS input file 



Figure A 1 . Flow chart of entire analysis procedure . 
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ABAQUS® input file 

beam.inp 


Figure A2. Flow chart of routine genbeam.f to generate finite element model. 
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program stress 
s S s start program 


INCLUDE ' aba_param. inc ' 

enable access to ABAQUS files 


r user specified 
input file 
stress.data 
contains loads, 
dimensions and 
. file names 


integer elnum, . . . 
real stress (6) . 
character . . . 

assign arrays 


read(nunit, *)word,bl, 

read data from user input file 


call gstart (kunit, output) 

open file for stress output 


LRUNIT ( 1 , 1 ) =8 
LRUNIT ( 2 , 1 ) =2 

CALL INITPF (FNAME , NRU, . 
CALL DBRNU(JUNIT) 

ABAQUS result file initialization 


stress output file 

beam.sig 


f ABAQUS result file > 
binary format 

^^^beam.fil^^^ 


CALL DBFILE (LOP, ARRAY, JRCD) 

read each record of the ABAQUS result file 


ELSE IF (KEY. EQ. 11) THEN 
call datin(head, 1 , . . . 

identify stress record and store data in array 


call sort (iunit, kunit, .. . 

sort stress data and write to output file 


call gend(iunit) 

close output file 


Figure A3. Flow chart of routine stress.f to retrieve stresses at the bottom of the beam. 

48 
















L: total specimen length 
s: span between supports 
s-L: span between inner load points 
t: specimen thickness 


(a) Two dimensional FE model of four point bending specimen with modeled rollers 
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Fig. A6a Comparison of Discrepancy between Beam Theory and 2D 
Plane-Stress non-linear FE analysis, S2/8552 4-Point Bend Configuration B2 
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Fig. A6b Comparison of Discrepancy between Beam Theory and 2D 
Plane-Stress non-linear FE analysis, 24-ply IM7/8552 4-Point Bend Config. B2 



B2, simple support 
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Fig. A6c Comparison of Discrepancy between Beam Theory and 2D 
Plane-Stress non-linear FE analysis, 36-ply IM7/8552 4-Point Bend Config. B2 
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